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Turbulent  mixing  of  multiphase  flow 

By  Y.-N.  Young 

J.  Ferziger,  F.  E.  Ham  and  M.  Herrmann 


1.  Motivation  and  objectives 

Mixing  of  multi-phase  fluids  is  common  in  diverse  research  fields,  and  understanding 
of  such  an  important  phenomenon  is  useful  to  a range  of  engineering  applications  such  as 
the  polymer  blender  or  control  and  designs  of  fluid  processing.  In  this  paper  we  adopt  the 
phase-field  modeling  approach  to  investigate  capillary  induced  effects  on  mixing.  Phase- 
field  modeling  has  been  applied  to  simulation  of  multi-phase  flow  (Chella  & Vinals  1996; 
Jasnow  & Vinals  1996;  Jacqmin  1999)  due  to  its  attractive  aspect  of  easy  numerical 
implementation.  More  recently  it  has  also  been  applied  to  simulations  of  solidification  of 
dendritic  alloys  (Zhao  et  al.  2003)  and  polymer  blenders  (Roths  et  al.  2002). 

On  macroscopic  scales,  the  interface  between  two  fluids  is  infinitely  sharp.  In  phase- 
field  modeling,  the  fluid  properties  are  modeled  (assumed)  to  change  smoothly  over  a 
small  layer  across  the  interface  on  mesoscopic  scales  (Chella  & Vinals  1996;  Jasnow  & 
Vinals  1996).  Such  a diffuse  interface  of  separation  is  first  considered  in  a liquid  near  or  at 
its  critical  point,  where  the  fluid  property  within  the  boundary  needs  to  be  modeled  and 
resolved  (Siggia  1979).  However,  phase-field  modeling  is  not  restricted  to  liquids  at  or  near 
the  critical  point  (Chella  & Vinals  1996;  Jasnow  & Vinals  1996;  Jacqmin  1999;  Young 
et  al.  2003).  In  fact,  one  can  certainly  envision  that  on  the  time  scale  of  macroscopic 
fluid  motions,  the  thermodynamic  potential  may  depend  on  the  gradient  of  the  order 
parameter,  which  designates  the  fluid  phase  (or  density).  Accordingly,  the  macroscopic 
jump  in  the  normal  stress  across  the  interface  due  to  surface  tension  can  be  modeled 
as  a mesoscopic  continuum  body  force  proportional  to  gradient  of  surface  potential. 
On  the  mescoscopic  scales,  such  a continuum  surface  tension  force  varies  smoothly  in 
the  transition  region  around  the  interface.  Naturally  such  a smooth  continuum  model 
of  fluid  interface  brings  great  numerical  convenience.  Recently,  phase-field  modeling  is 
made  exactly  energy-conserving  by  adding  a thermodynamical  constraint  (Jamet  et  al. 
2002).  We  did  not  include  this  constraint  in  this  work  as  the  correction  is  small  for  our 
cases  of  interest. 

An  interesting  finding  in  multiphase  fluids  stirred  by  chaotic  mixing  flow  is  the  self- 
similar drop  size  distribution  (Berthier  et  al.  2001;  Muzzio  et  al.  19916).  After  rescaling 
the  drop  size  distribution  with  respect  to  the  characteristic  size,  all  the  drop  size  distri- 
butions collapse  to  a universal  functional  form  for  various  combinations  of  parameters, 
such  as  viscosity,  capillary  number,  and  stirring  flow  strength.  Some  explanation  to  the 
self-similarity  is  provided  in  (Muzzio  et  al.  1991a, c),  and  it  is  not  obvious  if  similar 
collapse  of  size  distribution  can  also  be  found  in  a more  general  turbulent  flow,  where 
the  distribution  of  stretching  and  stretching  rate  is  randomly  distributed.  In  the  case  of 
stochastic  mixing  (Lacasta  et  al.  1995),  the  existence  of  a characteristic  domain  size  is 
an  indication  that  similar  collapse  of  drop  size  distribution  may  be  found.  However,  no 
report  on  the  self-similar  drop  size  distribution  is  given  in  Lacasta  et  al.  (1995).  Thus 
we  conduct  numerical  simulations  of  multiphase  fluids  stirred  by  two-dimensional  turbu- 
lence to  assess  the  possibility  of  self-similar  drop  size  distribution  in  turbulence.  In  our 


240 


Young  et  al. 

turbulence  simulations,  we  also  explore  the  non-diffusive  limit,  where  molecular  mobility 
for  the  interface  is  vanishing.  Special  care  is  needed  to  transport  the  non-diffusive  inter- 
face. Numerically,  we  use  the  particle  level  set  method  to  evolve  the  interface.  Instead  of 
using  the  usual  methods  to  calculate  the  surface  tension  force  from  the  level  set  function, 
we  reconstruct  the  interface  based  on  phase-field  modeling,  and  calculate  the  continuum 
surface  tension  forcing  from  the  reconstructed  interface. 

This  paper  is  organized  as  follows.  In  §2  we  formulate  the  coupled  Navier-Stokes-Cahn- 
Hilliard  system  and  present  results  from  direct  numerical  simulations  of  the  coupled 
system  stirred  by  turbulent  flows.  In  §3.1  we  formulate  the  multiphase  fluid  problem 
combining  the  particle  level  set  method  (for  interface  tracking)  with  phase-field  modeling 
(for  surface  tension  force).  Numerics  for  the  particle  level  set  method  are  summarized  in 
§3.2,  and  validation  is  provided  in  §3.3.  In  §3.4  we  present  results  of  drop  dynamics  from 
direct  numerical  simulations,  and  finally  we  provide  some  conclusion  in  §4. 


2.  Turbulent  mixing  in  the  Navier-Stokes-Cahn-Hilliard  system 

We  generalize  the  system  in  Berthier  et  al.  (2001)  to  include  the  effect  of  surface 
tension  on  stirring  flows.  The  phase-field  function  C denotes  the  phase  of  the  fluid: 
C = ±1  corresponds  to  fluid  phase  1 and  2,  respectively.  C is  described  by  the  Cahn- 
Hilliard  equation,  which  is  coupled  to  the  incompressible,  two-dimensional  Navier-Stokes 
equations  in  our  formulation. 

~ + u • VC  = /3MV2^,  (2.1) 

5 + u ■ Vu  = -VP  + i/V2u  - CY/ip  + F + Fd,  (2.2) 

at 

V • u = 0,  (2.3) 

where 

m = ^ F(C ) = \(C3  - l)2.  (2.4) 

In  equation  2.2,  F is  a random  forcing  at  wavenumber  k in  the  range  nx  < k/k0  < n2, 
where  ko  = 2n/l  and  l is  the  dimension  of  the  computation  domain.  Fd  = — Au  (A  = 0.1 
in  all  our  simulations)  is  the  drag  force  dissipating  energy  at  large  scales. 

The  positive,  constant  coefficient  M is  the  molecular  mobility,  u is  the  kinematic 
viscosity,  a and  /3  are  related  to  the  interfacial  properties:  The  interfacial  thickness  £ ~ 
t/oJP , and  the  surface  tension  coefficient  a = ^po^/aJS  = y/aj3cr0}  where  p0  is  the 
uniform  fluid  density  and  oo  is  the  surface  tension  force  when  \fa$  = 1.  We  can  non- 
dimensionalize  the  above  equations  with  respect  to  a characteristic  length  scale  l (of  the 
order  of  interfacial  thickness),  a characteristic  velocity  Vo,  and  the  corresponding  time 
l/v o-  The  dimensionless  equations  are 


— ■ + u • VC  = DV2  [C3  - C - V2C] , (2.5) 

^ + u • Vu  = -VP  + V2u  - ^CV  [C3  - C - V2C]  + F'  + F^,  (2.6) 

V • u = 0.  (2.7) 


The  dimensionless  coefficient  D = (3M/vqI  is  the  inverse  Peclet  number,  Re  = vqI/u  is 
the  Reynolds  number,  and  7 = \fafijv qv  is  the  capillary  number.  If  the  scaling  length 
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Figure  1.  Energy  spectrum  of  turbulent  velocity  (solid  line),  and  the  k s^3  inverse  cascade 
scaling  is  the  dashed  line,  ko  — ~jr  = 1/4. 


scale  l is  the  interfacial  thickness  (/  = £)  and  the  characteristic  velocity  is  the  “diffusion 
velocity”  across  the  interface  (u0  = then  the  Reynolds  number  Re  = (iM/v  and 

the  capillary  number  7 = a/f3Mv.  The  time  scale  for  diffusion  across  the  interface  is 
Tc  = £2 /PM.  The  characteristic  hydrodynamic  time  scale  is  ts  = 1/|S|,  where  |S|  is 
the  velocity  strain  rate.  It  is  convenient  to  control  the  interfacial  properties  using  both 
a and  /3,  thus  in  the  numerical  code  the  dimensional  equations  2. 1-2.4  are  implemented 
instead. 

For  the  following  simulation  results,  the  domain  size  of  the  double-periodic  computa- 
tion box  is  87 r x 8n,  and  the  numerical  resolution  is  51 22.  The  small-scale  forcing  range 
is  fixed  at  wavenumbers  5 < k/k0  < 15.  At  the  beginning  of  each  simulation,  an  initially 
circular  interface  of  radius  1.57T  is  placed  at  the  center  of  a two-dimensional  turbulent 
flow,  and  the  volume  average  of  the  concentration  is  ( C ) = 0.41  > 0.  The  energy  spec- 
trum of  the  turbulent  flow  is  the  inverse  cascade  with  a slope  of  —5/3  as  shown  in  figure  1. 
Throughout  the  simulations  the  energy  spectrum  is  little  affected  by  the  surface  tension 
force  from  the  interface. 

A numerical  challenge  in  turbulent  simulations  of  NS-CH  system  is  that  the  diffusion 
coefficient  has  to  be  small  to  avoid  phase  homogenization.  If  the  diffusion  coefficient  is 
too  large,  or  if  the  surface  tension  is  too  small,  drops  disappear  as  the  order  parameter 
C homogenizes  and  settles  to  the  phase  that  is  closest  to  the  mean  value.  Thus,  for 
a given  turbulent  strength  and  numerical  resolution,  we  optimize  the  number  of  drops 
by  minimizing  the  diffusion  with  enough  surface  tension  force  so  we  can  collect  drops  of 
various  sizes.  We  integrate  the  turbulent  flow  over  a long  duration  so  the  number  of  drops 
fluctuates  many  times  around  the  mean  value.  The  area  distributions  collected  from  three 
simulations  with  resolution  5122,  after  rescaled  with  respect  to  the  average  drop  area,  are 
shown  in  figure  2(a).  The  three  rescaled  area  distributions  collapse  reasonably  well  in  the 
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Figure  2.  Size  distribution  for  simulations  of  turbulent  NS-CH  flows.  Panel  (a):  Diamonds  are 
for  a = 2,  p — 200  and  M = 10~6.  Crosses  are  for  a = 2,  0 = 200  and  M = 5 x 10-7.  Asterisks 
are  for  a = 3,  = 300  and  At  = 10-6.  The  numerical  resolution  is  5122  and  u = 0.1  for  all  three 

simulations,  and  the  curve  is  ~ exp—  S/S*.  Panel  (b):  Diamonds  are  for  a — 2,  (3  = 200  and 
M — 2 x 10~7.  Asterisks  are  for  a = 1.5,  /3  = 150  and  M = 2 x 10-7.  The  numerical  resolution 
is  10242  and  v = 0.1  for  all  both  cases,  and  the  curve  is  ~ exp  —S/S*. 


range  0.3  < S/S*  < 20,  despite  the  poor  statistics  at  both  ends  of  the  distribution.  The 
diffusion  coefficients  (see  caption  in  figure  2)  that  we  used  in  the  numerical  simulations 
are  close  to  the  minimum  for  this  numerical  resolution  (5122).  10242  grid  points  are  used 
for  two  simulations  as  we  decrease  the  mobility  to  M = 2 x 10~7,  and  similar  scaling  in 
the  distributions  is  also  found  as  shown  in  figure  2(b).  Despite  the  limited  reliable  range 
in  both  figure  2,  the  general  trend  suggests  that  similar  scaling  may  also  exist  in  the 
tubulent  mixing  case. 


3.  Phase-field  modeling  with  zero  mobility  M = 0 

In  reality  the  mobility  M for  multiphase  flow  can  be  very  small,  and  the  equation 
for  the  phase  field  concentration  C reduces  to  the  advection  transport  equation.  This 
distinguished  limit  corresponds  to  the  non-diffusive,  immiscible  two-fluid  flow,  where  the 
interfacial  thickness  is  ideally  zero  and  thus  requires  special  numerical  treatment  in  the 
simulations. 

In  the  case  of  zero  diffusivity,  we  use  the  phase-field  modeling  as  a means  to  reconstruct 
the  interface  knowing  its  location.  We  use  particle  level  set  method  to  evolve  the  interface 
between  the  two  phases.  In  the  level  set  framework,  the  zero  of  the  scalar  level  set  function 
(j>  implies  the  location  of  the  interface.  To  retain  the  structure  of  the  interface,  we  assume 
that  the  Cahn-Hilliard  potential  ip  governs  the  interfacial  energetics,  and  reconstruct  the 
interface  for  a given  interfacial  thickness  £(=  y/a//3  as  in  §2).  In  §3.1  we  formulate  this 
system  in  terms  of  level  set  representation  of  the  interface  with  phase-field  modeling  for 
interfacial  structure  with  zero  molecular  mobility.  In  §3.2  we  briefly  discuss  the  numerics 
and  some  validation  is  given  in  §3.3.  We  then  present  results  from  direct  numerical 
simulations  of  multiphase  mixing  in  two-dimensional  turbulence  in  §3.4. 

3.1.  Formulation 

Basically  the  equations  are  exactly  the  same  as  those  in  §2  except  the  Cahn-Hilliard 
equation  is  now  reduced  to  an  advection  equation  for  C : 
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^ + (u-V)C  = 0.  (3.1) 

In  the  case  of  vanishing  molecular  mobility  across  the  interface,  the  advection  transport 
equation  3.1  for  concentration  C is  a numerical  challenge,  and  more  complication  arises 
in  the  evaluation  of  the  surface  tension  force.  Our  remedy  is  that,  despite  the  non- 
diffusiveness of  the  phase  field,  we  assume  that  the  interfacial  thickness  is  finite  based  on 
the  free  energy  if.  Naturally,  if  we  adopt  the  Cahn-Hilliard  surface  free  energy  for  if,  we 
expect  the  concentration  C to  behave  as  C ~ tanh(i//3/o|r  — ro|),  where  |r  — ro|  is  the 
distance  between  the  point  r and  the  interface  (located  at  ro).  A great  advantage  of  this 
approach  is  that  the  surface  tension  force  can  be  calculated  directly  from  the  phase-field 
modeling  as  a continuum  force  —CVif,  which  peaks  at  the  interface. 

Thus  instead  of  solving  equation  3.1  over  the  whole  computation  domain  for  C,  we 
solve  the  same  equation  for  a level  set  function  6 only  at  the  interface  (where  <j>  = 0),  and 
re-initialize  the  level  set  function  (f  to  a signed  distance  function  away  from  <j>  = 0.  We 
replace  equation  3.1  with  the  level  set  equation  valid  only  at  the  fluid  interface  <j>  — 0: 

dt<f>  + u-  V^  = 0|^=o.  (3.2) 

Away  from  <f  = 0,  inside  a band  of  width  ~ 6dx  around  cf  = 0,  we  re-initialize  <f  to  a 
signed  distance  function.  To  achieve  this  we  have  to  re-initialize  the  level  set  function  at 
each  time  step  by  solving  the  following  equation  to  a stationary  state 

fir* + «#»(*>)(!  V*|-1)  = 0,  (3.3) 

where  s gn((f0)  is  the  sign  of  0(x,  t = to)  defined  as 

S9n{<f>o)  = M + |v<M2d*2=’ 

for  a given  grid  spacing  dx.  The  normal  vector  and  the  mean  curvature  can  be  calculated 
in  terms  of  the  level  set  function  <f>\  n = V<£/|V<£|  and  k = V • n.  The  concentration 
phase-field  C is  then  reconstructed  from  the  level  set  function  (£,  as  previously  defined, 
is  the  interfacial  thickness): 

C = tanh(j)  if  \<f>\ < a£,  (3.5) 

c = ||i  if  M > a£, 

where  a is  a constant  in  the  range  2 < a < 5 so  that  the  concentration  C transitions 
smoothly  from  |C|  < 1 to  \C\  = 1 at  the  edge  of  the  band  around  the  interface.  Prom 
our  fully-resolved  simulations  of  drop  dynamics  presented  here,  we  conclude  that  the 
volume  integral  of  C remains  constant  at  all  time  provided  that  the  level  set  is  accurately 
capturing  the  interface. 

Numerically  there  are  various  ways  to  calculate  the  surface  tension  force  from  the  level 
set  function.  Physically  the  surface  tension  force  leads  to  a pressure  jump  across  the 
interface.  In  our  one-fluid  formulation  the  surface  tension  is  formulated  as  a body  force 
Fff  (=  —CVip  in  the  phase  field  modeling).  Ideally,  in  the  limit  of  zero  interfacial  thick- 
ness, the  surface  tension  force  is  non-zero  only  at  the  interface:  Fa  = (ann  + dsos)5((f), 
where  dsa  is  the  Marangoni  force,  with  s the  arclength  and  s the  unit  vectors  along  the 
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tangential  direction.  On  a discrete  numerical  grid,  however,  the  surface  tension  force  Fa 
has  to  be  smoothed  over  several  grid  spacings  to  avoid  un-controllable  numerical  oscil- 
lation. The  usual  cosine-delta  function  used  to  calculate  F„  on  a uniform  grid  (Sussman 
et  al.  1994)  leads  to  artificial  broadening  of  the  interface,  which  in  turn  leads  to  large 
unphysical  parasitic  currents  around  the  numerically  smoothed  interface.  In  our  phase 
field  modeling,  we  assume  an  interfacial  structure  based  on  the  the  Cahn-Hilliard  surface 
free  energy.  A continuum  surface  tension  force  is  then  derived  based  on  the  surface  free 
energy.  As  shown  in  (Jacqmin  1999),  this  continuum  surface  tension  force  gives  very  small 
parasitic  currents  compared  with  that  from  the  cosine-delta  construction  of  the  interface. 
Recently  a numerical  scheme  has  been  developed  to  minimize  the  parasitic  currents  (En- 
right et  al.  2002).  In  our  approach,  the  parasitic  currents  are  small  but  not  as  minimized 
as  in  Enright  et  al.  (2002).  The  turbulent  flow  in  our  resolved  simulations  of  drop  dy- 
namics is  much  larger  than  the  parasitic  currents,  and  thus  we  expect  the  constraint  in 
Enright  et  al.  (2002)  to  have  little  effects  on  the  dynamics  and  drop  area  distributions 
presented  in  §3.4.  The  formulation  we  proposed  for  reconstructing  the  interface  from  the 
phase-field  modeling  through  level  set  function  (equations  3. 2-3. 5)  is  somehow  similar 
to  the  splitting  of  particle  advection  from  particle  diffusion  in  the  random  walk  particle 
method  (Ghoniem  & Sherman  1985).  In  random  walk  particle  methods,  particles  can  be 
first  transported  by  the  velocity  field  at  their  positions,  and  then  perform  random  walks 
to  account  for  diffusion.  Such  a numerical  scheme  is  commonly  used  in  particle  tracking 
when  the  particle  diffusion  is  very  small.  Similarly  in  our  formulation  we  first  evolve 
the  interface  by  solving  the  advection  equation  only  near  the  interface  using  the  hybrid 
particle  level  set  method.  Once  the  location  of  the  interface  is  advanced  and  the  level  set 
is  re-initialized  as  a signed  distance  function,  we  construct  the  phase-field  C based  on 
the  steady  state  solution  to  equation  2.1  without  advecting  flow, 

C(C2  - 1)  - |V2C  = 0.  (3.6) 

Assuming  that  the  interfacial  structure  is  only  a function  of  the  signed  distance  to  the 
interface,  we  then  express  C in  terms  of  the  level  set  function  <p  in  equation  3.5. 

In  summary,  the  system  of  iso-density,  iso-viscosity  fluid  with  an  interface  separating 
the  two  phases  is  described  by  equations  2.2,  2.3,  and  3.2  with  the  surface  tension  force 
—CWip,  where  ip  is  the  Cahn-Hilliard  free  energy  and  the  concentration  C is  calculated 
from  equation  3.5. 


3.2.  Numerics 

The  level  set  equation  (3.2)  is  solved  using  the  hybrid  particle  level  set  method  (Enright 
et  al.  2002).  We  use  WENO  5th  order  scheme  to  calculate  the  flux  in  equation  3.2, 
and  the  TVD-3rd  order  RK  scheme  to  advance  the  level  set  in  time.  Following  Enright 
et  al.  (2002),  Lagrangian  particles  are  added  inside  a narrow  band  (rp)  around  the 
interface  and  are  transported  using  the  same  RK  3rd  order  scheme.  Before  and  after 
re-initialization  (equation  3.3),  particle  correction  (for  detail  see  Enright  et  al.  (2002))  is 
conducted  to  ensure  that  the  zero  level  is  not  shifted  numerically.  This  hybrid  level  set 
method  is  at  most  volume-conserving  as  long  as  particle  density  along  the  zero  level  is 
sufficient.  In  all  our  simulations,  we  use  at  least  42  particles  for  each  cell  within  the  (rp) 
band.  Re-seeding  of  particles  is  conducted  if  the  average  particle  density  is  found  to  be 
below  75%  of  the  initial  particle  density.  We  re-initialize  the  level  set  within  a second 
band  (rr)  enclosing  the  particle  band  (rp).  We  refer  readers  to  Enright  et  al.  (2002) 
for  detail  of  the  algorithm.  Standard  tests  have  been  performed  (Young  et  al.  2002)  and 
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results  are  identical  to  those  in  the  original  paper  (Enright  et  al.  2002).  To  couple  the 
level  set  equation(s)  to  the  fluid  solver,  we  advance  both  the  level  set  <t>  and  the  fluid 
velocity  u using  the  3-rd  order  explicit  scheme.  The  diffusion  part  in  the  momentum 
equation  is  incorporated  as  integrating  factors  in  spectral  space.  The  surface  tension 
force  is  calculated  from  the  interface  at  the  previous  time  step,  and  we  use  the  same 
velocity  fields  from  the  previous  time  step  to  advect  the  velocity,  interface,  and  particles. 

3.3.  Linear  capillary  waves 

To  demonstrate  that  the  surface  tension  force  is  properly  captured  in  the  phase-field  mod- 
eling, we  compute  the  oscillatory  frequency  of  the  capillary  waves  for  a two-dimensional 
drop.  The  initial  condition  for  (f>  is  the  signed  distance  function  perturbed  by  a sinusoid 
with  a small  amplitude  and  a wavenumber  n 

<j>(x,  t = 0)  = r — r0(l  + 0.02  sin(n<?)),  (3.7) 

where  ro  is  the  drop  radius,  r = |x-Xo|  is  the  distance  between  x and  the  drop  center  Xo, 
and  6 = tan_1((y— yo)/(x— xq)).  Analytically  (following  Rayleigh  (1892)),  the  oscillation 
frequency  of  capillary  waves  on  the  interface  separating  fluids  of  the  same  density  and 
zero  viscosity  (1/Re  = 0)  can  be  expressed  in  terms  of  surface  tension  coefficient  a,  drop 
radius  ro,  and  wavenumber  n of  the  perturbation 


We  obtain  the  linear  solutions  from  simulating  the  full  system  (described  in  section  3.1) 
with  small  non-linearity:  starting  from  the  initial  condition  u = 0 and  <j>  given  in  equation 
3.7,  the  non-linearity  remains  small  at  all  time  because  the  small  disturbance  is  damped 
by  viscosity.  In  all  our  simulations  in  this  section,  the  viscosity  u is  fixed  at  u = 0.1,  and 
the  perturbation  amplitude  decays  exponentially  at  rate  A with  a capillary  oscillation 
of  frequency  uj.  We  compute  the  capillary  frequency  u>  numerically  by  extracting  the 
oscillation  from  the  decaying  amplitude. 

We  first  calibrate  the  relationship  between  a and  (a,  (3).  Prom  §2  the  surface  tension 
er  is  proportional  to  With  the  ratio  (3/ a fixed  (the  interface  thickness  is  thus  fixed 
at  ~ 8 grid  spacings),  we  expect  the  surface  tension  a to  be  linearly  proportional  to 
a (or  •/?).  Table  3.3  shows  the  computed  oscillation  frequency  from  simulations  and  the 
corresponding  surface  tension  coefficient  o from  equation  3.8.  We  find  that  the  surface 
tension  a (computed  from  u>  in  our  case  via  equation  3.8)  indeed  varies  linearly  with 
y/afl. 

For  1 < n < 6 we  compute  the  oscillation  frequency  w from  the  simulations  with 
a = 0.9817  and  (3  = 89.1875  ( a = 94.32o-0  from  table  3.3).  As  shown  in  figure  3,  excellent 
agreement  is  found  between  the  computed  and  analytical  values  of  w for  all  1 < n < 6.  In 
addition  we  also  remark  that  the  parasitic  currents  are  as  small  as  10-6  with  a numerical 
resolution  1282  in  a square  box  of  size  2n  x 2ir. 

3.4.  Drop  dynamics  in  two-dimensional  turbulence 

Drop  dynamics  and  the  drop  break-up  in  simple  flow  configurations  has  been  extensively 
investigated  both  numerically  (Zaleski  et  al.  1995;  Cristini  et  al.  1998;  Tauber  et  al.  2002) 
and  experimentally  (Muzzio  et  al.  19916).  Results  from  experiments  (Muzzio  et  al.  19916) 
and  modeling  (Berthier  et  al.  2001)  on  passive  chaotic  mixing  of  immiscible  fluids  show 
that  the  statistics  of  drop  size  collapse  due  to  the  self-similar  breakup  and  coalescence 
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FIGURE  3.  Capillary  oscillation  frequency  w as  a function  of  n.  Symbols  are  simulation  data, 
and  line  is  from  analytical  result  for  a = 94.32cto  and  ro  = 7 r. 


Table  1.  Parameters  a and  /3  used  in  the  linear  simulations,  and  the  capillary  frequency  ui  com- 
puted from  the  time-varying  perturbation  amplitude,  n = 2,  v = 0.1,  and  numerical  resolution 
is  1282. 


a 

a 

U) 

cr/cro 

Case  1 

0.98 

89.18 

3.02 

94.32 

Case  2 

1.96 

178.37 

4.24 

188.64 

Case  3 

2.94 

267.56 

5.22 

282.96 

Case  4 

3.92 

356.75 

6.04 

377.28 

processes.  However,  there  are  only  a few  recent  direct  numerical  simulations  on  drop 
dynamics  in  turbulent  flows  (Tryggvason  et  al.  2001).  Due  to  the  difficulty  in  tracking  the 
interface  accurately  while  conserving  fluid  mass,  it  is  a great  numerical  task  to  capture  the 
statistics  of  drop  size  distribution  after  a series  of  breakups  and  coalescence  in  numerical 
simulations. 

Due  to  numerical  diffusion,  it  is  in  general  a great  numerical  challenge  to  conserve 
total  volume  for  each  fluid  using  level  set  methods.  The  hybrid  particle  level  set  method 
conserves  volume  almost  perfectly  in  standard  tests  (Enright  et  al.  2002),  and  is  the  only 
level  set  method  that  can  both  track  the  interface  accurately  and  minimize  the  mass 
loss  in  a controllable  fashion.  The  particle  level  set  method,  in  combination  with  the 
phase  field  modeling  of  the  surface  tension  force,  is  a good  tool  for  tracking  interface 
accurately  with  minimum  mass  loss.  Provided  that  enough  Lagrangian  particles  are  used 
in  tracking  the  interface,  the  mass  loss  can  be  reduced  to  a tolerable  level.  Thus  we  are 
able  to  accurately  evolve  the  dynamics  of  drops  in  a turbulent  flow,  and  our  main  goal  is 
to  examine  the  statistics  collected  from  DNS  dataset.  In  our  simulations,  less  than  1%  of 
maximum  volume  loss  is  found  for  all  cases  presented  here.  For  the  following  simulation 
results,  the  combination  of  a and  /?  is  chosen  so  that  the  interface  thickness  is  fixed 
(~  6 dx  at  least).  The  double-periodic  computation  domain  is  of  size  87T  x 871,  and  the 
numerical  resolution  is  5122. 
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t=0.6  1=1.9  t=13.1 


Figure  4.  Drop  deformation  and  break-up  in  two-dimensional  turbulence,  a = 2.95, 
ft  = 118.92,  v = 0.2  and  A = 0.1.  (a)  t = 0.6  (b)  t = 1.9  and  (c)  t = 13.1. 


We  fixed  the  range  of  small-scale  forcing  at  wavenumbers  5 < k/k$  < 15,  and  the 
large  scale  friction  F/  = -O.lu.  At  the  beginning  of  each  simulation,  an  initially  circular 
interface  of  radius  n is  placed  at  the  center  of  a two-dimensional  turbulent  flow.  The 
energy  spectrum  of  the  turbulent  flow  is  the  inverse  cascade  with  a slope  of  —5/3  as 
shown  in  figure  1.  Throughout  the  simulations  the  energy  spectrum  is  little  affected  by 
the  surface  tension  force  from  the  interface.  The  well-known  effect  of  bubbles/drops  on 
the  turbulent  energy  spectra  (Mazzitelli  et  al.  2003)  is  also  absent  here  because  there  is 
no  buoyancy /density  contrast  in  our  simulations. 

In  figure  4 we  demonstrate  how  the  initial  drop  gets  deformed,  stretched,  and  finally 
breaks  up  into  smaller  drops:  In  figure  4(a)  ( t — 0.6)  the  interface  is  highly  distorted,  and 
in  figure  4(b)  ( t = 1.9)  several  break-ups  lead  to  smaller  children  drops.  At  late  times 
(t  > 10)  the  drops  are  small  and  the  surface  tension  force  is  sufficient  to  overcome  the 
stretching  to  prevent  further  break-ups,  and  a statistical  equilibrium  is  reached  as  shown 
in  figure  4(c)  ( t = 13.1),  where  there  are  a lot  of  small  drops  of  irregular  shapes  and 
various  sizes. 

We  can  easily  calculate  the  total  length  of  the  interface  from  the  level  set  as 


where 


Length {$  = 0}  = / 

Jn 

-/ 

Jn 


\VH(<t>(x))\dxdy 

6(<f>(x))\y<t>\dxdy, 


* , i 


7r</>, 


- l + r + 5 Sin(-f) 


7 r 


(3.9) 


(3.10) 

(3.11) 


with  e = 3 ~ 4 dx  as  in  Sussman  et  al,  (1994).  Panel  (a)  in  figure  5 shows  the  evolution 
of  the  arclength  (scaled  to  the  initial  length),  and  panel  (b)  is  the  time-history  of  kinetic 
energy.  We  find  that  the  early  growth  in  the  length  is  accompanied  by  a similar  growth  of 
the  total  number  of  drops  as  several  break-up  processes  have  occurred  during  (0  < t < 3). 
After  t ~ 10  (about  2.5  eddy  turnover  times)  the  amplification  in  length  saturates  and 
fluctuates  throughout  the  simulations.  We  integrate  to  at  least  15  eddy  turnover  times 
to  attain  statistically  equilibrium  states.  The  interfacial  thickness  is  fixed  at  ~ 6.4 dx  for 
all  three  curves  in  figures  5.  With  the  same  small  scale  random  forcing  F,  large  scale 
dissipation  F^  and  viscous  dissipation  v = 0.2,  the  only  difference  between  the  three 
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Figure  5. 


(a)  Evolution  of  total  arclength  scaled  to  the  initial  value,  (b)  Corresponding 
evolution  of  number  of  drops. 


Table  2.  Averaged  number  of  drops  and  arclengths  for  the  three  simulations. 


a 

/? 

cr/ (To 

L/La 

N 

Case  1 

1.47 

59.45 

9.35 

7.51 

79.98 

Case  2 

2.94 

119.91 

18.71 

5.50 

39.69 

Case  3 

5.89 

239.83 

37.42 

3.66 

15.18 

simulations  is  the  surface  tension  coefficient:  a — 9.36er0,  <r  = 18.71cto,  and  a = 37.43<r0 
for  curves  1,  2,  and  3,  respectively.  (As  defined  earlier  cro  is  the  surface  tension  coefficient 
when  -y fafi  = 1.)  The  balance  between  stabilizing  surface  tension  force  and  stirring 
turbulent  flow  results  in  a critical  length  scale  (Lacasta  et  al.  1995) 


7T2  O 


c 2 u 


(3.12) 


where  u is  the  turbulence  characteristic  velocity.  Since  the  total  area  (volume)  of  each 
fluid  phase  has  to  conserve  in  two  (three)-dimensions,  we  can  estimate  the  dependence  of 
average  number  of  drops  (AT)  on  the  surface  tension.  For  a given  characteristic  turbulence 
velocity,  in  two-dimensions, 


NX 2 ~ total  area  of  fluid  phase  1 = constant  — > N rv  (J  (3.13) 

Similarly  in  three-dimensions, 

JVA3  ~ total  volume  of  fluid  phase  1 = constant  — ► N ~ cr-3//2.  (3.14) 

From  equations  3.13  and  3.14,  the  total  circumference  (area  or  arclength  in  two  and  three- 
dimensions,  respectively)  is  inversely  proportional  to  the  square  root  of  surface  tension: 
NXC  in  two-dimensions  and  N A2  in  three-dimensions  both  lead  to  ~ 1/ y/a.  These  scaling 
results  in  two-dimensions  are  consistent  with  the  time-averaged  values  of  our  simulation 
data  in  figures  5.  In  table  3.4  we  list  the  time-averaged  N and  total  arclength  for  all 
three  simulations. 

In  the  statistically  equilibrium  state  the  number  of  drops  fluctuates  around  the  mean 
value,  we  collect  enough  statistics  for  the  drop  area  distribution  by  integrating  over  a 
long  period  of  time.  Following  Berthier  et  al.  (2001)  and  Muzzio  et  al.  (19916),  we  rescale 
the  drop  area  distribution  by  the  mean  area  (calculate  by  taking  moments  of  the  area 
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FIGURE  6.  Size  distribution  for  the  three  cases.  Diamonds  are  for  case  1,  asterisks  are  for  case  2 
and  crosses  are  for  case  3.  Open  circles  are  for  a different  forcing  function  but  same  parameters 
as  in  case  1.  The  dashed  line  is  ~ exp  —(S/5*)  . 

distribution  functions  as  in  Muzzio  et  al.  (19916)).  Since  the  effective  diffusion  coefficient 
is  assumed  to  be  zero  in  our  turbulence  simulations,  we  cannot  calculate  the  mean  area 
based  on  the  scaling  of  the  characteristic  domain  size  L*  and  D*  as  in  Berthier  et  al. 
(2001).  The  rescaled  area  distributions  collapse  as  shown  in  figure  6.  One  additional 
simulation  data  set  (open  circles)  with  a different  interfacial  thickness  and  turbulent 
random  forcing  is  added  in  figure  6.  For  this  additional  simulation,  the  surface  tension 
a = 18.71ero  and  the  interfacial  thickness  £ = 5.0dx  and  the  turbulent  random  forcing  is 
at  wavenumbers  5 < k/k0  < 7.5,  narrower  than  the  other  three  simulations.  However,  it 
seems  that  the  collapse  works  reasonably  well  regardless  of  the  differences  between  the 
simulations. 
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4.  Conclusion 

We  have  extended  the  phase-field  modeling  to  the  non-diffusive  case.  Surface  free 
energy  from  the  phase-field  modeling  is  used  to  calculate  the  surface  tension  force  based 
on  the  reconstructed  interface  using  the  level  set  function.  Numerically  the  non-diffusive, 
infinitely  sharp  interface  has  to  be  smoothed  over  a few  grid  points.  Some  interfacial 
structure  is  always  assumed  in  the  numerics,  and  the  resultant  surface  tension  force 
inherits  the  modeled  structures.  One  popular  way  to  reconstruct  interface  in  level  set 
method  is  the  cosine-delta  function.  However,  no  physical  meaning  can  be  related  to 
such  a construction,  and  as  a result,  the  parasitic  currents  are  undesirably  large.  We 
replace  the  cosine-delta  function  with  the  smooth  interfacial  structure  from  the  stationary 
solution  in  phase-field  modeling,  and  we  calculate  the  surface  tension  force  from  the 
corresponding  surface  free  energy.  The  only  assumption  we  make  in  this  approach  is 
that  the  reconstructed  interface  is  a function  of  the  signed  distance  to  the  interface.  We 
assume  that  “some”  molecular  dynamics  acts  on  such  a fast  time  scale  that  the  interfacial 
thickness  remains  the  same,  even  when  the  flow  is  vigorously  acting  on  and  around 
the  interface.  The  collapse  of  drop  area  distribution  in  our  two-dimensional  turbulence 


250  Young  et  al. 

simulation  confirms  our  conjecture  that  the  self-similar  drop  area  distribution  exists  for 
more  general  velocity  fields. 

The  particle  level  set  method,  combined  with  phase-field  modeling,  has  great  promise 
for  large  scale,  three-dimensional  numerical  simulations  of  multiphase  flow.  One  of  our 
future  goals  is  to  generalize  our  simulations  of  two-dimensional  drop  dynamics  to  three 
dimensions.  We  also  plan  to  include  density  and  viscosity  contrast  for  multi-fluid  using 
the  phase-field  modeling.  Finally,  we  are  also  investigating  how  the  phase-field  modeling 
can  be  applied  in  simulating  micro-fluidic  multiphase  flows  carried  out  in  experiments 
(Truesdell  et  al.  2003). 
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